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1.  Introduction 


There  are  many  models  of  ballistics  trajectories.  The  high-resolution  6-degree-of-freedom 
(6-DOF)  models  require  many  computations  and  a  small  time  increment.  The  modified  point 
mass  models  ignore  the  spinning  of  the  round  to  reduce  the  computational  requirements. 

Selecting  a  model  to  use  for  ballistics  estimation  or  tracking  requires  tradeoffs  between  system 
accuracy  and  computation  expense.  For  example,  using  a  6-DOF  model  for  a  tactical  (real-time) 
system  is  not  currently  possible;  however,  for  experimental  work,  it  is  possible  to  interpret  the 
data  using  a  6-DOF  model.  The  most  basic  type  of  polynomial  filter  is  a  straight-line  predictor. 
These  are  accurate  over  short  distances,  but  they  will  ignore  any  curvilinear  behavior  (typically 
needed  for  intercept  prediction),  and  their  perfonnance  can  suffer  when  the  prediction  time  is 
increased.  In  a  similar  fashion,  higher-order  polynomial  predictors  will  diverge  as  prediction 
time  increases.  Reducing  the  fidelity  of  the  model  saves  time  and  expense  associated  with 
computational  requirements,  but  does  so  by  sacrificing  accuracy.  The  purpose  of  this  report  is  to 
present  ballistics  models  (simpler  than  the  6-DOF  model)  that  can  be  used  to  model  the 
dynamics  in  an  extended  Kalman  filter  (EKF).  The  EKF  linearizes  the  dynamics  at  the  system 
operating  point  and  then  proceeds  as  a  Kalman  filter  (KF). 

Intercept  systems  require  the  tracking  of  the  ballistics  threat  and  the  interceptor.  When  there  is  a 
limited  amount  of  control  authority  available,  estimation  must  be  very  accurate.  As  control 
authority  increases  and  terminal  guidance  sensors  improve,  it  is  possible  to  use  models  of  lower 
fidelity  to  estimate  the  trajectory  of  a  round  (target  or  interceptor).  Other  situations  also  make 
reduced-order  modeling  possible.  For  the  nonlinear  effects  that  are  small,  it  may  be  possible  to 
remove  them  from  the  dynamic  model.  This  includes  estimation  with  a  high  sampling  rate  where 
the  time  between  observations  is  small  and  trajectories  where  the  drag  coefficient  does  not 
change  too  much. 

A  trajectory  can  be  adequately  modeled  by  a  differential  equation.  Using  the  initial  conditions, 
an  expected  projectile  path  can  be  generated.  Range,  velocity,  position,  and  direction 
measurements  can  be  used  to  improve  the  estimate  of  the  projectile’s  path.  The  differential 
equation  and  measurement  process  are  combined  to  fonn  a  KF.  The  differential  equation  models 
the  physics  of  the  trajectory,  while  the  measurement  is  used  to  update  the  parameters  of  the 
equation  through  least-squares  estimation. 
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2.  One-Dimensional  (1-D)  Case 


First  it  is  helpful  to  consider  the  1-D  case.  The  1-D  case  is  useful  because  it  forms  the  basis  of 
understanding  three-dimensional  (3-D)  trajectories;  for  some  direct-fire  situations,  it  provides  an 
excellent  model  for  projectile  analysis.  In  situations  where  time  until  impact  is  of  interest  for 
high-velocity  rounds,  this  model  can  provide  timing  infonnation. 

2.1  Basic  Equations 

If  a  force,  F,  is  acting  on  a  body  and  the  resistance  to  the  force  is  proportional  to  the  velocity 
squared,  we  have  the  following  straightforward  differential  equation  from  Newton’s  law: 

m  —  —  F—bv 1 .  (1) 

dt 

2  F 

If  k  =  — ,  then  the  equation  can  be  written  as  follows: 

£  =  -V-*’).  (2) 

dt  m 

Solving  this  equation  involves  a  separation  of  variables,  as  shown  in  the  following: 


1  ,  b 

- -  dv  = - dt  . 


Z  /  z 

v  -  k 


Using  a  partial  fractions  expansion, 


P  |  g 


(v2  -  k2)  v-  k  v  +  k 


Watching  your  p’s  and  q’s  leads  to  the  following  differential  equation: 

1  ,  1  1  .  ,  b  , 

— ( - )dv  = - dt .  ( 

2k  v—k  v+k  m 

Integrating  this  results  in  the  following: 

— (ln(v  -  k)  -  ln(v  +  k))  - — —t  +  c.  ( 

2k  m 

Taking  the  exponential  of  each  side  and  using  properties  of  exponents  results  in  the  following: 


From  this,  the  velocity  can  be  found  as  a  function  of  time,  as  shown  in  the  following: 
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(8) 


v(t)  —  k 


1  +  ce~a/ 

1  -  ce~at  ' 


If  a 


2 kb 

m 


and  c  = 


v0~k 

v0  +  k 


where  v0  is  the  initial  velocity,  the  velocity  will  asymptotically 


approach  k.  This  equation  can  be  used  to  find  the  terminal  velocity  of  a  dropped  object  by  letting 
F  be  the  force  due  to  gravity.  A  large  value  of  F  could  be  used  to  simulate  the  launch  of  the 
projectile.  After  launch,  the  value  of  F  would  be  0  unless  a  rocket  or  missile  is  being  modeled. 


In  the  case  of  flat  trajectory  over  a  short  time,  the  only  force  acting  on  the  object  is  drag.  Solving 
this  same  equation  for  F  =  0  yields  the  following  result: 

-«-rr-  <9) 

- 1 - t 

v0  m 


This  expression  asymptotically  approaches  0.  Assuming  v(t)  has  been  measured,  it  is  possible  to 
use  the  previous  equation  to  find  the  initial  velocity  using  the  following: 


1 

vo=“r~r 

v(t)  m 


(10) 


Given  knowledge  of  two  velocities,  v(h)  and  v(i2) ,  it  is  also  possible  to  calculate  the  value  of 
drag  ( b ).  In  the  following,  note  that  b  includes  the  effects  of  air  pressure  (everything  but  speed) 
and  t  is  the  difference  between  the  two  times: 


b  =  ( 


1 

v(t2) 


1  m 
v(q)  t 


(ID 


Position  can  be  found  by  adding  distance  traveled  to  the  original  position.  The  distance  traveled 
is  found  by  integration  of  velocity  and  is  described  by  the  following  expression: 


m 

~b 


ln( 


1  b 

- H - 

v0  m 


o- 


m 

~b 


ln(— ). 


(12) 


The  distance  traveled  is  the  natural  logarithm  of  a  linear  function.  Using  this  expression,  an 
upper  bound  for  the  lateral  distance  traveled  can  be  established.  It  is  also  possible  to 
approximate  the  initial  velocity  needed  to  attain  a  specified  distance.  This  equation  provides 
another  constraint  the  analyst  can  use  in  fitting  data.  The  major  problems  with  this  formulation 
are  the  1-D  restriction  and  the  assumption  of  constant  drag.  Typically,  drag  is  a  function  of 
Mach  number.  Modeling  of  drag  has  resulted  in  universal  drag  curves.  These  curves  give  the 
overall  shape  of  drag  as  a  function  of  Mach  number.  For  individual  rounds,  these  curves  are 
multiplied  by  a  number,  typically  called  a  form  factor,  to  adjust  the  universal  curve  to  adequately 
fit  the  particular  round.  Drag  changes  dramatically  in  the  region  of  Mach  1;  thus,  nonlinear 
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behavior  can  be  a  concern  in  this  region.  The  form  factor  also  allows  the  use  of  drag  curves  of 
the  same  projectiles  to  be  moved  up  or  down  for  the  round  of  interest.  This  adjustment  accounts 
for  slight  perturbations  of  the  projectile’s  shape  and  mass.  The  drag  curves  for  similar  shells  are 
assumed  to  be  similar;  the  form  factor  allows  this  to  be  incorporated  into  the  model  of  motion. 
While  the  incorporation  of  a  more  complex  drag  model  would  increase  the  fidelity  of  the  model, 
it  would  preclude  a  closed-form  solution.  Incorporation  of  more  complicated  drag  models  will 
require  the  equations  to  be  solved  numerically.  The  closed-form  model  allows  a  computation  of 
time  until  impact  and  other  useful  quantities  associated  with  time  or  distance  along  the  trajectory. 
The  selection  and  design  of  a  1-D  model  will  depend  on  the  application. 

2.2  Measurements 

An  EKF  can  be  designed  for  the  1-D  case.  For  information  on  Kalman  filtering,  see  Gelb  (/)  or 
Maybeck  (2).  The  equations  for  an  EKF  using  the  notation  given  by  Gelb  follow. 

System  nonlinear  dynamics  plus  plant  noise  q~N(0,Q): 


X  =  f(x(t),t)  +  <7(0 . 

(13) 

The  observation  equation  v~N(0,R): 

Zk  =  K  (x(tk ))  +  vk  . 

(14) 

Initial  conditions,  normal  distribution: 

x(0)~A(x(0),P(0)). 

(15) 

The  covariance  propagation: 

Pk  =  F{x{t),t)Pk_ i  +  Pk_xF{x{t),t)  +  Q(t ) . 

(16) 

The  gain  due  to  an  observation: 

Kk  =  PkHk  (xk  )[ffk  (x,  )PkH'k  (x, )  +  Rk  ]_I . 

(17) 

Change  in  the  state  due  to  observation: 

K  =xk  +  Kk(zk  -hk(xk)) 

(18) 

Updated  state  covariance  via  observation: 

(19) 

Linearized  time  step: 

F(x(t),t)  =  ^  ^  x(0  =  xk  . 

dx(t) 

(20) 
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Relinearized  observation: 


H{xk) 


SJh  ( x(h )) 

dx(tk) 


\x{t)  =  Xk  . 


(21) 


The  state  is  a  set  of  parameters  that  allows  the  differential  equation  to  be  solved.  First,  define  the 
state  vector,  X,  as  the  position  (X  i),  the  speed  (X2),  and  the  drag  coefficient  (X  3).  Given  this 

set  of  parameters,  it  is  possible  to  propagate  the  trajectory  using  the  following: 

Xx(k)  =  Xi(k-l')  +  X2(k-V>dt  +  \x22  (k-^Xidt2  ’  (22) 


x^^x^-^+x^-Ox/*  ’ 


(23) 


and 


X3W  =  X3 


(24) 


where  k  indicates  the  time  step.  Assuming  dt  is  small,  dt 2  will  be  close  to  0.  These  equations 
can  be  written  as  follows: 


(° 1 

. 0  ) 

f  \ 

X, 

x(k)=x(k-l)+ 

0  0  x/k-1) 

x2 

0  0 

0  J 

yXi) 

(25) 


The  3x3  matrix  is  referred  to  as  F  and  captures  the  change  as  a  function  of  time.  Assume 
P(k—  1)  is  the  state  covariance  and  using  E  to  represent  expectation,  the  new  state  covariance  is  as 
follows: 

E(X  (k)x'(k))  =  (I  +  F)E(x  (k  - l)X(k  - 1))(7  +  F)  (2fi) 

*  P(k  - 1)  +  FP(k  - 1)  +  P(k  - 1  )F, 

since  the  expectation  of  %  (k-  l)j^  (k  - 1)  is  the  state  covariance  P(k—  1).  This  formula  allows  a 

closed-form  means  of  predicting  the  variation  of  a  future  state  from  the  current  conditions.  By 
using  equations  22-26,  is  possible  to  predict  the  striking  velocity  and  also  have  statistical  bounds 
on  the  variation  of  the  striking  velocity.  Covariance  propagation  is  an  alternative  to  Monte  Carlo 
simulation  of  the  system.  The  change  of  the  system  and  the  covariance  of  the  system  in  time  are 
referred  to  as  state  propagation  and  state  covariance  propagation.  When  an  observation  of  a  state 
variable  or  combination  of  state  variables  occurs,  the  state  can  be  updated  in  a  least-squares 
manner. 
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Observations  will  change  the  perception  of  the  state.  The  current  example  will  be  developed  to 
include  the  effect  of  location  or  distance  measurement  and  velocity  or,  in  this  case,  speed 
measurement,  on  both  the  state  and  the  state  covariance.  The  measurement  error  of  the  location 
and  speed  are  assumed  to  be  known  to  a  reasonable  degree  of  accuracy. 


Typically,  the  time  between  observations  is  fixed  and  dt  is  constant;  however,  this  need  not  be 
the  case.  For  asynchronous  measurements,  the  value  of  dt  can  be  adjusted  to  meet  the  situation. 
If  the  interval  between  measurements  becomes  large,  it  is  prudent  to  update  the  covariance 
matrix  between  observations.  When  a  single  measurement  is  made,  the  state  is  updated  via 
recursive  least  squares.  The  measurement  is  expressed  in  terms  of  the  state  variables  and  will  be 
signified  with  the  symbol,  h(X(t )) ;  the  observation  matrix  will  be  the  partial  of  this  with  respect 

to  the  state  and  represented  by  H.  In  terms  of  the  state,  a  position  observation  is 
(1  0  0)  X;  thus,  the  //matrix  is  (1  0  0).  For  a  velocity  measurement,  the  observation 

corresponds  to  (0  1  0)  X.  In  both  these  cases,  the  measurement  is  a  linear  function  of  the  state; 
when  the  measurement  is  a  nonlinear  function,  the  partial  of  measurement  is  used  to  define  H  at 
each  time  step.  First,  the  gain  matrix,  k,  will  be  defined  based  on  elements  of  the  state 
covariance  matrix,  P,  and  the  measurement  error,  as  shown  in  the  following: 


(  Pn) 


k  = 


Pn 

VPJ 


/VCT 


(27) 


In  this  situation,  i=  1  for  a  position  measurement,  and  i  =  2  for  a  velocity  measurement,  a 
represents  the  standard  deviation  of  the  measurement,  and  the  subscript  will  indicate  a  position 
or  velocity  measurement.  The  distance  variance  will  differ  from  the  velocity  variance  and  may 
change  from  measurement  to  measurement.  The  state  update  due  to  the  observation  is  as 
follows: 

x(+)  =  x{-)  +  k{z-Xi).  ^8) 

In  this  equation,  z  is  the  measurement,  and  the  -  sign  indicates  before  the  observation  update. 
Notice  that  the  terms  of  the  k  matrix  are  directly  proportional  to  the  corresponding  terms  of  the 
covariance  matrix.  The  observational  update  always  reduces  the  state  covariance,  as  shown  in 
the  following: 

P(+)  =  P(-)-D—^,  (29) 

P„  +CT, 


and 


P.P.,  . 
ij  ik 


(30) 
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Together,  propagation  and  observation  updates  form  a  KF.  Since  a  KF  is  recursive,  it  needs  to 
be  initialized.  Typically,  the  experimenter  initializes  the  filter;  it  is  important  to  have  a  good 
initial  estimate  of  the  state  and  the  state  covariance. 

The  drag  coefficient  is  more  realistically  represented  as  a  function  of  Mach  number.  As  the 
projectile  changes  speed,  its  drag  coefficient  changes.  It  is  usually  infeasible  to  generate  drag 
curves  for  every  projectile;  a  characteristic  curve  is  found  for  a  particular  projectile.  Given  a 
similar  projectile,  the  curve  is  moved  up  or  down  by  incorporating  a  form  factor  into  the 
equation.  In  equation  25,  X,  would  become  the  fonn  factor.  The  drag  coefficient  would  be 

included  in  the  F  matrix  as  a  factor  of  the  square  of  the  speed.  Note  that  state  propagation  is  a 
numeric  solution  of  the  differential  equation.  In  the  case  of  an  extended  filter,  the  dynamics  are 
linearized  at  each  time  step.  The  1-D  model  can  be  used  when  the  projectile’s  motion  stays  on  a 
line  or  possibly  when  timing  is  the  important  information.  This  model  can  also  be  used  to  find 
the  distance  along  a  known  trajectory.  Note  also  that  the  1-D  model  can  be  enhanced  by 
including  the  effects  of  gravity  to  form  a  simple  two-dimensional  (2-D)  model. 


3.  Two-Dimensional  Case 


The  2-D  case  is  more  complex  because  the  velocity  term  is  an  interaction  of  both  the  height  and 
range  terms  (cross  range  is  ignored).  In  a  sense,  this  interaction  steals  velocity  from  the  range 
dimension  as  gravity  causes  the  vertical  velocity  to  asymptotically  approach  its  terminal  velocity. 
In  time,  the  direction  of  motion  will  align  itself  with  gravity.  Only  drag  and  gravity  are 
considered;  the  equations  are  as  follows: 

x  =  -cos(0)  —  v2  ,  (31) 

m 

y  =  g~  —  v2  sin(0)  ,  (32) 

m 

v2  =  x2  +  y2  ,  (33) 


and 


0  =  atan(y  /  x)  . 


(34) 


To  solve  this  system,  the  initial  conditions  must  be  stated.  The  previous  model  (equations  3 1— 
34)  has  been  used  to  model  submunitions  being  released  from  a  cargo  round.  Assume  the 
submunitions  are  expelled  in  the  range  dimension  so  that y  -0-0  and  x  =  v0 .  Choosing  the 

expulsion  velocity  allows  these  equations  to  be  solved  numerically.  The  previous  equation  was 
realized  in  SIMULINK  and  solved  numerically  therein.  For  horizontally  launched  submunitions, 
the  lateral  speed  approaches  0  as  the  terminal  velocity  is  attained.  The  value  b  is  based  on  the 
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drag  coefficient;  a  discussion  of  this  can  be  found  in  U.S.  Army  Special  Text  9-153  (3)  or 
Sabersky  et  al.  ( 4 ). 

The  value  of  drag  at  velocities  less  than  Mach  1  is  fairly  constant.  If  the  velocity  crosses  Mach 
1,  a  more  complicated  model  for  drag  should  be  devised.  A  2-D  model  for  prediction  can  be 
useful  for  endgame  guidance  applications  where  the  time  to  hit  is  low  and  the  relative  position 
information  from  the  sensor  is  good.  The  previous  model  (equations  3 1-34)  has  been  used  to 
model  the  dispersion  of  submunitions  (5).  Equations  3 1-34  can  be  used,  but  it  would  be  more 
accurate  to  reduce  the  order  of  the  equations  discussed  in  the  next  section. 


4.  Three-Dimensional  Case 


Flight  dynamics  are  most  accurately  represented  by  6-DOF  models.  These  models  are  nonlinear 
and  must  be  solved  numerically.  There  are  no  closed-form  solutions,  and  the  numerical  solutions 
require  many  computations  at  each  iteration;  it  is  difficult  to  develop  a  real-time  filter  based  on  a 
6-DOF  model.  Many  estimators  for  flight  dynamics  use  simpler  models  focused  on  the 
parameters  of  interest.  For  example,  if  position  is  of  interest,  then  it  is  possible  to  ignore  some  of 
the  computations  associated  with  attitude  and  reduce  the  complexity  of  the  model. 

Point  mass  models  and  modified  point  mass  models  offer  a  simplification  of  the  6-DOF  models 
that  provide  excellent  position  accuracy.  In  three  dimensions,  a  drift  term  must  be  added  to  the 
model;  this  term  captures  a  projectile’s  motion  orthogonal  to  range  and  altitude.  Drift  is  caused 
by  the  interaction  of  spin  and  yawing  motion.  Spin  can  also  be  modeled,  and  its  decrease  is  in 
proportion  to  the  current  spin  rate  and  the  speed  of  the  round.  Obtaining  precise  knowledge  of 
aerodynamic  coefficients  can  be  difficult,  and  even  with  this  knowledge,  there  can  be  round-to- 
round  variation.  It  is  prudent  to  include  a  form  factor  to  capture  the  variation  of  aerodynamic 
coefficients. 

For  a  linear  model,  the  concept  of  state  is  used  to  find  a  Markov  representation  of  the  system. 

The  state  is  a  vector  that  incorporates  the  information  needed  to  propagate  forward  to  the  next 
time  of  interest.  Nonlinear  dynamics  are  often  modeled  by  linearization  of  the  model  followed 
by  choosing  a  time  step  that  does  not  result  in  nonlinearities  of  a  problematic  magnitude.  A  14- 
dimensional  state  model  will  be  discussed  (this  follows  the  presentation  from  excerpts  from  a 
portion  of  an  unidentified  report).  The  state  is  shown  in  appendix  A.  It  is  possible  to  reduce  this 
state  model  to  eight  (position,  velocity,  drag,  and  lift)  or  seven  (position,  velocity,  and  drag) 
dimensions.  Also,  universal  curves  will  be  used  to  model  aerodynamics  associated  with  drag, 
lift,  and  spin.  It  is  assumed  that  a  user  with  more  knowledge  of  a  particular  round  can  replace 
the  modeled  dynamics  with  higher  fidelity  models,  but  these  will  represent  the  default  models. 
The  14  dimensions  include  three  for  position,  three  for  velocity,  two  for  wind  speed,  and  one 
each  for  muzzle  velocity,  azimuth  bearing,  elevation  bearing,  drag  constant,  lift  (or  drift) 
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constant,  and  speed  of  sound  at  sea  level.  A  seven-dimensional  (7-D)  model  can  be  formed  by 
using  position,  velocity,  and  drag.  An  eight-dimensional  model  would  add  lift  or  perhaps  muzzle 
velocity  to  the  previous  state.  The  number  of  elements  in  the  state  can  be  chosen  based  on  the 
purpose  of  the  model. 

The  dynamics  can  be  thought  of  in  tenns  of  the  acceleration  of  the  round.  Additive  components 
are  grouped  as  drag,  drift,  gravity,  Magnus  effect,  and  Coriolis  terms.  A  description  of  modified 
point  mass  models  can  be  found  in  STANAG  4355  ( 6 ).  As  a  frame  of  reference,  the  North-Up- 
East  system  attached  to  the  projectile  launch  point  on  the  surface  of  the  earth  will  be  used. 

The  drag  term  is  a  function  of  air  pressure,  velocity,  round  diameter,  and  aerodynamic 
coefficients.  For  the  calculation  of  drag,  velocity  is  squared;  thus,  velocity  is  the  most  important 
tenn  in  the  drag  equation.  The  drag  coefficient  is  calculated  from  a  fourth-degree  polynomial  of 
Mach  number,  with  coefficients  based  on  universal  drag  values.  There  are  other  methods  that 
can  be  used  to  model  the  drag  coefficient.  The  dynamics  associated  with  drag  expressed  in  tenns 
of  the  state  variables  are  as  follows: 


In  this  case,  Cd  is  found  using  the  universal  drag  curve.  V  represents  the  speed  of  the  projectile 
relative  to  the  ground.  The  symbol  xu  is  the  fonn  factor  or  the  factor  that  allows  the  universal 

drag  value  to  be  moved  up  or  down.  The  variable  A  represents  the  air  pressure.  Air  pressure 
changes  with  height,  so  A  =  f(x2).  The  standard  atmosphere  model  is  typically  used  to  find  air 

pressure,  although  meteorological  data  can  be  used  if  available.  Additional  factors  relating  to  the 
physical  characteristics  of  the  round  can  also  be  included. 

The  lift  term  in  modified  point  mass  models  is  most  accurately  described  using  the  yaw  of 
repose.  The  yaw  of  repose  represents  the  projectile’s  average  yaw.  For  this  example,  the  yaw  of 
repose  will  not  be  modeled,  and  the  lift  tenn  will  be  orthogonal  to  the  projectile  velocity  and  the 
gravity  vector.  (See  STANAG  4355  for  a  method  to  estimate  yaw  of  repose.)  The  term 
associated  with  the  lift  will  be  as  follows: 


In  this  formula,  C,  is  the  lift  term.  In  the  present  situation,  it  will  be  calculated  from  the 
universal  lift  curve.  The  form  factor  xu  is  used  to  adjust  the  lift  curve  for  the  current 

application.  V  is  the  speed  of  the  projectile.  The  universal  lift  curve  was  developed  using  rounds 
with  high  spin  rates;  other  methods  for  rounds  with  a  low  spin  rate  perhaps  can  be  ignored. 
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The  Magnus  force  will  be  ignored  in  the  dynamics  used.  This  force  is  important  for  predicting 
impact  time.  With  regular  observations  of  projectile’s  position,  this  term  will  not  have  a  large 
impact  between  observations.  The  omission  of  the  Magnus  force  term  will  increase  the  model 
uncertainty.  This  tenn  was  ignored  in  the  modified  point  mass  model  discussed  by  Bradley  (7). 


Gravity  needs  to  be  included  in  the  model  dynamics.  The  force  of  gravity  is  directed  to  the 
center  of  the  earth.  The  magnitude  of  gravity  changes  over  the  earth  as  a  function  of  latitude. 
The  following  formula  is  an  approximation: 


g  =  go(1_£i  cos(2 L)), 

(37) 

g0  =  9.80665, 

(38) 

and 


Si  = 


.0026, 


(39) 


where  L  is  the  latitude  at  the  point  of  launch.  In  terms  of  the  state  variable,  the  gravity  vector  is 
as  follows: 


g  =  g 


xJR 

I-O^+jc2)'5 1(2  R) 


x3/R 


(40) 


where  R  is  the  radius  of  the  earth. 


Coriolis  force  is  a  factor  when  using  an  earth-fixed  coordinate  system.  The  rotation  of  the  earth 
is  7.2921  x  10  ^  =  Q  .  Let  Q r  =  Geos (Latitude)  and  Qv  =  Qsin  (Latitude),  then  the  Coriolis 

effect  can  be  written  as  follows: 


x4 

Qv*6 

*5 

=  -2 

-nrx6 

nxx5-nyx4 

(41) 


The  dynamics  discussed  are  summarized  in  appendix  B.  These  can  be  used  as  the  basis  of  an 
EKF.  If  a  7-D  state  is  being  used,  the  dynamics  associated  with  lift  can  be  omitted. 


In  situations  where  guidance  is  required,  a  model  of  the  spin  is  necessary.  STANAG  4355 
proposes  the  following  differential  equation: 


P= 


npd^vcs 

K 


(42) 
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where  p  is  spin,  p  represents  the  air  pressure,  d  is  the  diameter  of  the  round,  cs  is  the  spin  drag 
coefficient,  and  Ix  is  the  moment  about  the  spin  axis.  C  is  a  function  of  Mach  number.  This 

equation  could  be  used  as  the  dynamics  of  a  separate  EKF  to  estimate  spin,  assuming  there  are 
some  measurements  of  spin  available.  In  conjunction  with  the  previous  position  model,  the 
combination  of  the  two  models  for  position  and  spin  could  estimate  much  of  the  information 
desired  about  a  projectile’s  flight.  The  combination  of  these  two  models  would  not  capture 
yawing  motion.  This  may  not  be  an  issue  if  the  yaw  is  small,  thus  allowing  the  system  to  be 
modeled  using  a  simplified  or  reduced  state  estimation. 

The  dynamics  can  be  used  to  develop  a  state  estimator.  The  theory  is  discussed  in  Gelb  (7).  The 
development  of  a  KF  varies  based  on  the  dynamics  being  modeled;  thus,  it  is  possible  to  have 
many  different  KF  estimating  the  same  quantities.  Differences  in  KF  are  due  to  the  state 
propagation  equations  used.  For  nonlinear  problems,  an  EKF  is  a  good  first  choice.  An  EKF  is 
based  on  the  same  theory  as  a  KF  but  uses  a  linearized  version  of  the  nonlinear  dynamics  at  each 
time  step.  If  the  time  step  is  made  too  large,  the  EKF  may  not  be  appropriate,  and  the  model 
may  need  to  include  more  terms  related  to  the  nonlinear  dynamics.  Another  alternative  is  using 
particle  filters.  Particle  filters  do  not  require  the  propagation  of  the  state  uncertainty;  this  benefit 
is  offset  by  the  need  to  propagate  a  number  of  candidate  states  forward  in  time  and  then  calculate 
uncertainty  based  on  the  distribution  of  candidate  states.  In  addition  to  nonlinear  dynamics 
associated  with  the  state  transition,  there  can  be  nonlinearities  associated  with  the  observations. 
An  observation  needs  to  be  expressed  in  terms  of  the  state  variables.  If  position  is  part  of  the 
state  and  the  distance  from  an  object  to  the  projectile  is  observed,  then  the  observation  is  a 
nonlinear  function  of  the  state.  Nonlinearities  in  the  state  dynamics  and  nonlinearities  in  a 
measurement  expressed  in  terms  of  the  state  variables  are  mitigated  through  using  an  EKF  by 
linearization  around  the  current  value  of  the  state. 

To  develop  an  EKF,  proceed  with  the  following  steps.  First,  decide  on  the  dynamics  to  capture 
by  the  EKF.  Next,  the  dynamics  need  to  be  put  in  the  form  x  -Fx  +  Gu,  where  x  is  the  state  of 

the  system  and  F  represents  the  change  in  x  over  a  time  step.  The  variable  u  represents  inputs  to 
the  system,  and  the  matrix  G  describes  how  these  inputs  affect  the  system.  Another  use  for  the 
Gu  term  is  to  introduce  uncertainty  associated  with  unmodeled  dynamics.  After  these  tasks  have 
been  completed,  the  process  must  be  initialized.  Then,  as  observations  become  available,  these 
are  incorporated  into  the  filter  through  least-squares  estimation.  As  time  progresses,  the  state 
and  its  covariance  will  propagate  forward  in  time.  The  state  propagation  equations  have  been 
discussed.  Only  position  and  velocity  change.  The  equations  for  an  EKF  using  the  notation 
given  by  Gelb  can  be  found  in  section  2.2. 

Using  the  previously  discussed  dynamics  and  the  EKF  equations,  an  estimator  for  the  trajectory 
of  a  projectile  can  be  designed.  Note  that  the  matrix  to  be  inverted  has  the  size  of  the 
observation  covariance;  typically,  this  is  smaller  than  the  state  covariance.  The  one  issue  not 
discussed  is  the  observation  equation.  A  GPS  sensor  gives  position  so  the  observation  matrix 
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consists  of  a  3  x  N  matrix  of  zeros,  with  ones  in  the  positions  corresponding  to  location,  so  that 
Hxk  gives  the  position  in  terms  of  the  state  at  the  Mi  observation. 


Many  observations  are  nonlinear;  radar  typically  gives  the  range  to  the  target.  Assuming  the 
radar  is  located  at  (0,0,0)  and  the  first  three  values  of  the  state  correspond  to  position,  the  radar 
observation  in  terms  of  the  state  is  r  =  (x2  +  xi  +  x32) 5 .  If  range  was  the  only  available 


observation,  the  linearized  observation  matrix  would  contain  one  row,  the  first  three  columns 


would  be  H(x)  = 


and  the  other  entries  would  be  zeros.  If  azimuth  and  elevation 


angles  are  given,  they  can  be  expressed  in  terms  of  the  position,  and  then  the  proper  partial 
derivatives  can  be  found  and  included  as  extra  rows  of  H. 


Typically,  most  of  the  effort  associated  with  an  EKF  goes  into  finding  the  matrix,  F(x(t),t) ,  the 
partial  of  the  system  dynamics  with  respect  to  the  state.  This  matrix  can  be  complex  and  is 
typically  recomputed  each  time  step  or  observation.  If  the  F  matrix  does  not  change  too  quickly, 
it  is  possible  to  process  several  time  steps  before  recomputing  the  matrix.  Also,  if  some  of  the 
partials  are  numerically  small,  it  is  reasonable  to  drop  these  terms  from  the  matrix.  It  is  possible 
to  have  many  different  formulations  of  an  EKF  for  the  same  problem.  The  interplay  between 
desired  accuracy  and  computational  speed  detennines  the  final  fonn  of  the  EKF. 


5.  Example  for  the  3-D  Case 


In  this  section,  an  EKF  with  a  7-D  state  using  GPS  measurements  will  be  discussed.  The  state 
variables  will  consist  of  position,  velocity,  and  ballistics  coefficient.  The  dynamics  for  lift  and 
Magnus  effect  will  be  ignored.  The  state  variables  used  from  appendix  A  will  be  1-6  (position 
and  velocity)  and  12  (ballistics  coefficient).  Also,  wind  effects  will  be  assumed  to  be  zero  to 
simplify  the  simulation.  A  GPS  measurement  gives  position;  and  since  position  is  part  of  the 
state,  this  results  in  an  H  matrix  with  a  3-D  identity  matrix  followed  by  zeros  in  positions 
corresponding  to  the  fourth  through  seventh  state  variables.  The/ equation  represents  the 
dynamics  as  follows: 


fr* 4< 

(43) 

fl=X5  ’ 

(44) 

ii 

(45) 

d  Vx4  -  ~~  ~  2G  vx6  , 

(46) 

12 


fs  =  -x7AkdVx5  -  g(l  -  (Wj)  )  +  2Q-VX6  > 


2  R, 


gx  3 


fe  =~x7AkdVx 6  -  —  -2Qxx5  +2Qvx4, 


(47) 

(48) 


and 


fi=  o. 


(49) 


Using  equations  43-49,  the  F  matrix  is  found  by  taking  the  partials  with  respect  to  each  state 
variable  (see  equation  20).  Recall  that  both  kd  and  V  are  functions  of  x4,  jc5,  and  x6.  In  equations 

43-49,  V  is  the  speed  of  the  projectile  and  the  drag  coefficient,  kd  ,  is  a  function  of  Mach 
number,  which  is  a  function  of  speed.  Also  note  that  air  pressure,  A,  is  a  function  of  altitude,  x2 . 
Assume  the  F  matrix  starts  out  as  a  7  x  7  matrix  filled  with  zeros.  The  following  identifies  the 
nonzero  elements: 


734=1, 


^25=1, 


736=1: 


F 


41 


-g 

R. 


(50) 

(51) 

(52) 

(53) 


F42  =-x7kdVx4 


dA 

dx0 


-^44  —  XqA 


VxA^-  +  kdxA  —  +  k,V 


dxA 


8x 


A45  —  x7Ax4 


VdAL  +  k  BV 


4  J 

\ 


V  dx5 


dx 


5 ) 


(54) 

(55) 

(56) 


F46=~x7Ax4 


vaAL+kA 1 


V  8X6 


dx. 


-2Q,, 


6j 


7^51  = 


F 47  =  -AkdVx4  , 
g  xi 


2  \.5 


2Re  (Xx  +  X3  ) 

p  -  -x  k  Vx  — 

r52  ~  xlKdvx5  ~ 

OXn , 


(57) 

(58) 

(59) 

(60) 
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F  =  — 

/53 


g  *3 

2 Re  (• * 2  +*32)'5  ’ 


654  —  X1  ^X5 


{v^  +  k,F 

{  dx4  dx4  j 


F55  —  x7A 


V 


^56  ~  XjAx5 


Vx5d^L  +  kdX5j-  +  k 
OX  c  OW 

+  2Q 


8k,  ,  5F 3 


y—+kd  — 

dx6  dx 


6  J 


F57  —  AkdVx5 

P  8A 

F62  =  ~x7kdVx6  — 


dx2  ’ 


77  - _ £ 

63  “  R  ’ 


6*64  —  X?Ax6 


7*65  —  X7Ax6 


dk. 

a 

+  kd - 

dx4 

&4 

dkd 

+  kd - 

dx5 

dx  5, 

+  2  Q 


-2Q 


‘X  9 


7*66  ~  X-jA 


Vx6^L  +  kdx6^  +  kdV 

dx,  dx. 


and 

F67  =  "^^6  ' 

The  following  infonnation  is  also  needed  to  obtain  numerical  values  for  the  F  matrix: 

V  =  (x4  +x 52  +Xg)'5  , 

5x,.  (X4  +  x52  +x62  )'5  *  ’  ’ 


a0  =1.223 , 
ax  -  1.071e  -  4 . 


J  =  a0e  01X2 


(61) 

(62) 

(63) 

(64) 

(65) 

(66) 

(67) 

(68) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 

(75) 

(76) 
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-01*2 


(77) 


and 


8A 

—  —  —  aGaxe 

dx2 

g  =  9.80665(1  -  .0026  cos(2  Lot)) ,  (78) 


Re  =6356766.  (79) 


The  drag  coefficient  is  calculated  using  a  fourth-degree  polynomial  of  Mach  number.  Mach 
number  is  the  speed  divided  by  the  speed  of  sound.  Notice  that  in  this  formulation,  the  partial  of 
the  speed  of  sound  with  respect  to  height  is  not  included  in  the  F  matrix. 


s0=  340.3  temperature  =  59°  F , 

(80) 

xsl  =height(launchabovesealevel ) , 

(81) 

cv  =  2.26e-5 , 

(82) 

j  =  j0(1-cv(jc,,  +v2))-5, 

(83) 

4 

kd  =  XC',J7'  > 

i= 0 

(84) 

dkd  ( .  ,-_i  'j  xi 

8Xi  It?  Jsr’ 

(85) 

Q  =  7.2921e  -5  , 

(86) 

Qr  =Qcos  (/at), 

(87) 

Qv  -  Qsin(/at) . 

(88) 

In  this  example,  the  measurements  to  be  used  are  assumed  to  be  estimates  of  position  from  a 
GPS  receiver,  so  the  observation  matrix  H  is  expressed  as  follows: 


H  = 


1  0 
0  1 
0  0 


0  0  0  0 
0  0  0  0 
10  0  0 


0 

0 

0 


(89) 


When  multiplied  by  the  state  vector,  this  matrix  will  select  the  three  elements  associated  with 
position  as  the  three  observations.  For  these  observations,  there  are  no  issues  with  nonlinearities. 
The  next  step  is  to  use  the  preceding  equations  to  track  a  projectile.  The  number  of  observations 
per  second  can  be  varied  as  a  design  parameter.  A  6-DOF  model  will  be  used  to  define  a 
trajectory. 


15 


The  initialization  of  the  estimator  will  be  discussed.  The  initial  value  of  the  state  can  be  set 
using  the  launch  position,  launch  speed,  azimuth,  and  elevation  of  the  launch.  The 
R  matrix  is  the  covariance  of  the  observation.  For  the  present  case,  this  can  be  found  from  the 
specifications  of  a  GPS  receiver.  In  other  situations,  this  could  involve  building  a  model  of  the 
sensors  and  then  performing  a  sensitivity  analysis  to  get  an  accurate  estimate  of  the  R  matrix.  In 
many  cases,  the  R  matrix  will  change  due  to  changes  in  the  target  sensor  geometry.  One 
R  matrix  will  be  used  for  this  EKF.  The  Q  matrix  will  be  used  to  model  the  shortcomings  of  the 
dynamic  model  used  in  the  EKF.  Typically,  a  guess  is  made  for  the  Q  matrix.  This  guess  is  then 
adjusted  through  the  information  gained  by  repeated  adjustments.  The  process  of  tuning  an  EKF 
or  KF  is  finding  a  reasonable  Q  matrix. 

The  EKF  can  be  used  to  approximate  system  performance  measures  based  on  state  covariance. 
The  state  covariance  matrix  P  is  available  at  each  time  step  and  can  be  used  to  derive  measures 
of  effectiveness.  This  type  of  analysis  can  be  used  to  determine  the  observation  rate  and  sensor 
quality  needed  to  meet  performance  specifications.  This  covariance  analysis  does  not  require  the 
extensive  use  of  replications  required  by  Monte  Carlo  simulations  and  is  closer  to  a  closed-form 
solution.  As  an  example,  consider  an  observation  rate  of  10/s.  Assume  the  observation 
covariance  is  a  diagonal  matrix,  with  4,  9,  4  on  the  diagonal;  also,  let  the  Q  matrix  be  diagonal, 
with  ones  in  positions  one  through  six  and  .0001  in  position  seven.  For  a  mortar  round  shot 
north  at  49°  elevation,  an  aggregate  measure  of  the  state  uncertainty  is  the  trace.  Figure  1  shows 
the  trace  of  the  state  covariance  matrix  as  a  function  of  time  in  units  of  0.01  s. 
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Figure  1.  State  covariance  trace  vs.  time. 
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From  this  graph,  it  is  seen  that  at  ~5  s  the  steady  state  behavior  is  achieved.  The  band  is  the 
result  of  covariance  stochastically  expanding  due  to  both  the  Q  matrix  and  the  F  matrix  until  an 
observation  is  made  and  then  instantly  diminishing  as  a  result  of  the  new  information.  This  can 
be  seen  in  the  blowup  of  figure  1  presented  as  figure  2. 


Figure  2.  Expansion  of  figure  1 . 

To  determine  how  well  this  filter  works,  data  from  a  6-DOF  model  was  obtained,  and  the 
previously  described  7-D  EKF  was  used  to  track  the  data.  In  figure  3,  the  blue  is  the  6-DOF 
track,  and  the  red  is  the  output  of  the  EKF.  It  is  difficult  to  distinguish  the  two  curves. 

The  reduced  dynamic  model,  coupled  with  the  observations,  tracks  the  6-DOF  data  with  limited 
error. 


6.  Conclusions 


It  is  realistic  to  simplify  the  model  when  the  nonlinearities  being  ignored  do  not  adversely  affect 
the  estimation  or  result  in  errors  that  are  within  tolerances.  Also,  when  observations  are 
available,  these  observations,  due  to  their  accuracy,  may  correct  the  estimator  to  the  extent  that 
including  some  of  the  dynamics  is  not  worth  the  computational  expense;  that  is,  with  a  high  data 
rate,  it  is  often  possible  to  use  a  simple  model  and  achieve  acceptable  results.  In  linear  systems 
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Figure  3.  True  vs.  estimated  trajectory. 

theory,  the  effects  of  unmodeled  dynamics  are  referred  to  as  model  mismatch  error.  There  are 
uncertainty  terms  that  can  account  for  the  unmodeled  dynamics;  this  leads  to  an  increase  in  the 
covariance  associated  with  the  state  variables.  It  is  not  possible  to  present  all  possible  reduced 
state  models;  the  intent  was  to  show  how  it  can  be  done  and  discuss  some  of  the  issues  associated 
with  designing  an  EKF. 

It  is  possible  to  break  the  problem  into  a  position,  velocity  model,  and  spin  model.  In  many 
cases,  using  these  two  separate  models  will  provide  an  acceptable  alternative  to  using  6-DOF 
models  for  target  tracking. 

The  accuracy  of  a  ballistics  estimator  needs  to  be  considered  in  system  design.  There  is  a 
tradeoff  between  accuracy  and  computation  requirements.  For  a  kinetic  energy  round,  a  1-D 
model  that  incorporates  gravity  may  suffice.  The  observations  available  to  the  system  will 
influence  the  choice  of  dynamics.  If  there  are  many  high-quality  observations,  a  simple  linear 
predictor  may  suffice.  In  this  situation,  the  nonlinear  effects  may  be  small  between  observations. 
The  matrix  F  representing  the  partials  of  the  state  dynamics  with  respect  to  the  state  variables  is 
computationally  the  most  difficult  term.  Terms  in  the  F  matrix  with  relatively  small  magnitudes 
are  to  be  considered  as  candidates  for  omission.  For  target  interceptor  systems,  the  overall 
system  performance  is  determined  by  the  accuracy  of  each  estimator  (target  and  interceptor)  in 
conjunction  with  the  control  authority. 
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Appendix  A.  The  State  for  a  Three-Dimensional  System 


Location 


*1 

Position  north 

*2 

Position  up 

x3 

Position  east 

Velocity 

*4 

Speed  north 

*5 

Speed  up 

*6 

Speed  east 

Wind  Velocity 

x7 

Speed  north 

xs 

Speed  east 

Xg 

Muzzle  speed 

Vo 

Azimuth  angle  to  target 

*n 

Elevation  angle  to  target 

xl2 

Drag  constant 

*13 

Lift  constant 

*14 

Speed  of  sound  at  sea  level 
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Appendix  B.  Change  in  the  State  Variables 

A=X 4‘ 

fl  =X5- 

A  =X6- 

A  =  ~XnAkdV (X4  -  V  )  -  <X  -  *8  )  -  m yX6  ■ 

v  Ke 

fs=-xnAidVxs -gil~(x'  +X;)S)  +2ClA  . 

1Ke 

f6=-xi2Akdv(x6-xs)-^L(x4-xi)-^L-2nxx5+myx4 . 

V  Kc 

/v-/l4=0. 
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